Conformational dependence of a protein kinase phosphate transfer reaction 
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Atomic motions and energetics for a phosphate transfer reaction catalyzed by the cAMP- 
dependent protein kinase (PKA) are calculated by plane-wave density functional theory, starting 
from structures of proteins crystallized in both the reactant conformation (RC) and the transition- 
state conformation (TC). In the TC, we calculate that the reactants and products are nearly isoen- 
ergetic with a 0.2 eV barrier; while phosphate transfer is unfavorable by over 1.2 eV in the RC, with 
an even higher barrier. With the protein in the TC, the motions involved in reaction are small, 
with only P^ and the catalytic proton moving more than 0.5 A. Examination of the structures 
reveals that in the RC the active site cleft is not completely closed and there is insufhcient space 
for the phosphorylated serine residue in the product state. Together, these observations imply that 
the phosphate transfer reaction occurs rapidly and reversibly in a particular conformation of the 
protein, and that the reaction can be gated by changes of a few tenths of an A in the catalytic site. 

PACS numbers: 



INTRODUCTION 

Protein kinases regulate many biological processes by 
transferring a phosphate group from adenosine triphos- 
phate (ATP) to the sidechains of particular serine, threo- 
nine, or tyrosine residues. The bulky, charged phosphate 
group alters the conformation and function of the tar- 
get protein [l|, |2|. Different kinases recognize different 
primary sequence motifs surrounding the residue to be 
phosphorylated, in a highly regulated fashion [^|j,|5|,y{. 
Structural studies have revealed several conformational 
changes, such as closing of the active-site cleft, the pack- 
ing of the activation loop, and rotation of the C-helix, 
which arc often implicated in controlling the activity of 
protein kinases J2|. The reasons for such control are clear, 
but no answer has been provided to such questions as: 
"How closed is closed?" , or "is this particular conforma- 
tion of the activation loop 'good enough' for phospho- 
rylation to occur?" Quantum chemistry is required to 
objectively answer these questions. 

The extent of conformational heterogeneity in a cova- 
lent protein reaction was first quantified in a series of 
experiments monitoring the temperature-dependent re- 
binding of CO to myoglobin after flash photolysis [2|. 
Agmon and Hopfield created a concise phenomenologi- 
cal model describing this situation, using transition state 
theory to describe the vibrational reaction, and a dif- 
fusive coordinate which describes the protein conforma- 
tion and modulates the reaction barrier to the vibrational 
transition |3, Q . Moving beyond the phenomenological 
model requires specifying both the conformational het- 
erogeneity and the sensitivity of the reaction barrier to 
this heterogeneity. Careful structural analysis [l3 and 
quantum chemistry calculations [ll| showed that the re- 
action barrier heterogeniety is indeed a reasonable conse- 
quence of observed structural heterogeniety at the myo- 
globin active site. A reevaluation of a wide variety of 



myoglobin data also shows that a distinction between 
diffusive, solvent-controlled conformational motions and 
Arrhenius transitions which are independent of the sol- 
vent dynamics is well supported by experiments [l^l . 

In this work, we explore the conformational sensitivity 
of the protein kinase reaction in two experimentally de- 
termined structures of PKA, one crystallized with ATP 
and the protein kinase inhibitor (PKI), which we refer 
to as the reactant conformation of PKA, or RC [13j. 
The other conformation is obtained by crystalizing with 
a transition-state analogue, a non-reactive ADP-AIF3 
and a mimic of PKI which is both shorter and has a 
phosphate-accepting serine instead of an inert alanine at 
the reactive position [lj| . We refer to the protein confor- 
mation in this case as the transition state conformation 
of PKA, or TC. Although differences in conformation 
between the two structures are ~ 0.5 A, we find qual- 
itatively different energetics of reaction, which lead us 
to conclude that conformational motion of the protein 
kinase is rate-limiting to the overall phosphate transfer 
reaction. 



METHODS 

Initial equilibrium geometries are generated from the 
protein data bank entries lATP, (RC) and 1L3R, (TC), 
both crystal structures of PKA. Initial guesses for the 
complete reactant and product structures were obtained 
by homology-modeling the terminal phosphate of ATP 
(P^) and side chains of the serine and catalytic aspartic 
acid (D-166) into appropriate positions. All atoms which 
were modeled in this way were allowed to move in subse- 
quent geometry optimization steps. Four different model 
sizes, containing 82, 88, 244 and 263 atoms, were con- 
structed. The 82-atom minimal system has been defined 
by Valiev et al. |ia| , and includes the Mg2-tri phosphate 



and its immediate interaction partners, G-52, S-53, K- 
72, D-166, K-168, N-171, D-184, and the substrate serine 
residue. Two water molecules from the crystal struc- 
ture are included to complete the six-fold coordination 
of Mg in the 88-atom system. The 244-atom system con- 
tains a significantly larger shell around the reaction cen- 
ter and more of the ATP molecule (including part or 
all of residues G-52, S-53, F-54, K-72, D-166, L-167, K- 
168, P-169, E-170, N-171, D-184, F-185, G-186, and the 
substrate backbone from residue 18 to 22). The atomic 
coordinates of backbone carbon and nitrogen atoms for 
each system were taken directly from the experimental 
crystal structures. Non-backbone bond lengths and most 
angles are derived from amino acid templates to facilitate 
comparison of energies between structures. Protons were 
added where needed. Residues were truncated by replac- 
ing carbon atoms with protons and adjusting the new 
C-H bond length accordingly. Full atomic coordinates 
are provided in the supporting material. 

Atomic interactions are described by the VASP de nsity 
functional theory (DFT) and a plane waves basis |l6l . 
Ultrasoft pseudopotentials of the Vanderbilt form 
and a PW91 generalized gradient approximation func- 
tional are used. A 270 eV plane wave cutoff, appropri- 
ate for the pseudopotentials, is applied. This cutoff was 
increased to 300 eV to verify insensitivity of results to 
choice of basis set. The reaction pathways arc computed 
with periodic images separated by at least an 8 A vac- 
uum layer. A periodic box size of 19 x 21 x 16 A'^ is used 
for the 82- and 88-atom clusters, 20 x 24 x 21 A^ for the 
244-atom cluster and 27 x 24 x 23 A^ for the 263 atoms 
structure. 

Reactant and product geometries are calculated by op- 
timizing the geometry of the clusters with a conjugate 
gradients algorithm. Saddle points connecting these sta- 
ble geometries were found using the nudged elastic band 
(NEB) method [l9j. In this method, images are gener- 
ated by linear interpolation between the optimized reac- 
tant and product structures, and DFT is used to calcu- 
late forces on each atom of each image. The images are 
then geometry optimized subject to harmonic forces be- 
tween the images which force them to be equally spaced 
along a minimum energy pathway between rcactants and 
products. The climbing image modification to the NEB 
method [20, |2l| was used so that the highest energy im- 
age along the band converges directly to the saddle point, 
thus increasing the accuracy of the energy barrier with 
fewer images. Refinements to the barriers for small sys- 
tem changes were computed with the dimer saddle point 
finding method |22|. 

The 82-atom system contains the same 59 uncon- 
strained atoms as Valiev et al., while the waters added 
to make the 88-atom system were allowed to move. The 
244-atom system has only 42 moving atoms (the gamma 
phosphate, Mgl, Mg2, their coordinated water, and the 
sidcchains of S-21, D-166, and K-168). 



In order to test the sensitivity of the 244-atom result 
to system size changes and details of the electrostatic 
boundary conditions and periodic box size, we added the 
adenosine ring to make a 263 atom system, and observed 
changes to the barrier for TC of less than 15 meV. 



RESULTS 
conformational dependence 

Fig. ^ shows the energetics of the 244-atom systems. 
The reaction is endothermic in RC by 1.2 eV, while in 
TC, the reaction is nearly isocnergetic with a barrier of 
only 0.20 eV. The reaction pathway of TC, shown in Fig. 
^, shows only small motions during the reaction; only 
two atoms move more than 0.5 A, and the catalytic base 
moves only 0.07 A. Apparently, in the TC structure, the 
atoms around the active site are correctly arranged for 
both the reactants and products state. In contrast, 1.2 
eV endothermicity in RC indicates sub-optimal geome- 
tries in the product state. 

At the transition state, the P7O3 has a planar configu- 
ration, halfway between the ADP and the serine O7, and 
proton transfer from the serine to the Asp- 166 occurs af- 
ter the phosphate transfer. Transferring the proton to 
an oxygen atom on the P-y phosphate was found to be 
unfavorable by more than 0.5 eV. These observations are 
in agreement with previous DFT studies [l5j. 

One possibile spurious origin for the low barrier in TC 
is that both the reactants and products arc destabilized 
because the active site has collapsed around the AIF3, 
just as we expect RC to be energetically biased towards 
the reactants. Two observations force us to the conclu- 
sion that this is not occuring. 

Examination of the reaction pathway geometries re- 
veals two clear reasons for the endothermicity in the TC 
structure. First, the reaction cleft in RC is approximately 
1 A more open in RC than in TC, so the protein is un- 
able to maintain the full octahedral coordination of both 
magnesium ions in the products state. Fig. |21 shows this 
distance between the N-171 oxygen and the closest P/3 
oxygen to be 5.1 A in the RC conformer. In the TC 
structure, the reaction center is compressed, reducing 
this distance to 4.2 A. Second, it appears that there is 
more space for the phosphorylated serine residue in TC 
than RC. 

Fig.^ maintains the relative energies of RC and TC, 
and shows that TC is only 160 meV higher than RC. The 
errors contributing to this difference are likely several 
tenths of an eV, but it is gratifying that the reactants 
energies are so similar in the two conformations. 
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FIG. 1: (a) Energy barrier for the 244-atoni system in the RC 
and TC conformers. With the additional constraining atoms, 
the RC reaction becomes very unfavorable, while the TC con- 
formation phosphorylates. (b) Reaction path for phosphate 
transfer in PKA for a 244-atom TC structure. The additional 
constraining atoms on the outside of the cluster, as compared 
to the 82-atom system in Fig. EI a-) i reduce spurious motion 
during the reaction. For example, the central Mg^+ ion re- 
mains in place with the proper coordination throughout this 
reaction. 



geometry constraints 

Another question arises from the relatively small set 
of moving atoms allowed in the calculation shown in Fig. 
^ One might expect the geometry-optimized reaction 
pathway to be independent of initial conformation in the 
limit of the entire protein being allowed to move. Several 
observations force us to the conclusion that this is not 
correct. 

First; energetic minimization of both RC and TC to- 
wards the reactants state, using a classical molecular dy- 
namics potential (which should be perfectly valid for the 
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TABLE I: The distances in A that particular atoms moved in 
the four reaction pathways. The distance from the reactants 
minimum to the transition state is indicated in parenthesis 
for the 244-atom TC structure. 



reactants state) causes only one third of the difference in 
the distance shown in Fig.|21to dissappear, even with no 
solvent present and every atom in the protein allowed to 
move during the optimization. 

Second, Table I, shows much smaller motions are 
needed for reaction in TC than in RC. The Mg++ ions, 
coordinated water molecules, catalytic base, and lysine 
all move - 0.1 A in TC, and 0.5 to 1.0 A in RC. Even 
the serine residue, which gains a phosphate group, moves 
only a quater A during the course of reaction. A re- 
lated observation is that the change in force on the frozen 
atoms during reaction is about half as large for TC than 
RC (not shown). 

To directly check the effect of constraint on the results, 
and guided by knowlege of the forces on frozen atoms, we 
relieved the constraint on the beta phosphate of ATP. 
The relative energy of RC and TC went up by 350 meV, 
while the maximum change in force dropped by a factor 
of ten, and was uniformily distributed across a dozen 
boundary atoms. The exothermicities of RC and TC 
changed by less than 0.1 eV, while optimal geometries 
differed by ~ 0.1 A, suggesting that the templates used to 
build the Mg-triphosphate coordination sphere differed 
from the quantum mechanical potential by ~ 0.1 A in 
several places. 

The conclusion that TC provides a structure with 
happy reactants as well as the potential to create happy 
products with minimal motions is robust. 



dependence on system size 

Valiev et al. [15j reported a barrier of 0.5 eV for this 
reaction, when using either a GGA or B3LYP functional 
and local basis set on an 82-atom version of the RC struc- 
ture, with only a pair of atoms fixed at the boundary 
of each molecular fragment. This result is intermedi- 
ate between our RC and TC results, and could indicate 




FIG. 2: The 82-atoin system in the TC structure has a lower 
barrier than the RC structure and a favorable product state 
because the reaction center is compressed. In the product 
and transition states, the central Mg atom can remain octa- 
hedrally coordinated in the TC conformer. In the RC con- 
former the reaction cleft is opened and the same Mg atom is 
forced to break one bond, making the reaction unfavorable. 
The distance between the coordinating oxygen atoms on the 
Asnl71 and ADP groups, which is a measure of reaction cen- 
ter size, is increased from 4.3 A in the TC structure to 5.0 A 
in the RC structure. 



computational details provide variability as large as the 
conformational differences which we cite. 

To check this posibility, we duplicated their 82-atom 
system and choice of constraints, and show the results of 
our calculation in Fig. 13 for both RC and TC. Consis- 
tent with their result, wc find a 0.5 cV barrier and an 
exothermicity of 0.2 eV in the RC, indicating a robust- 
ness of the result to the various truncation proceedures, 
basis sets, periodic boundary conditions, and functionals 
used the calculations. Unfortunately, the result is also 
insensitive to the starting conformation of the protein, in 
contrast to the results on the 244-atoni systems. 

Table I shows that the motions of the atoms are much 
larger than in the 244-atom calculation, and compari- 
sion of Figs. ^ and|2l3 shows a qualitatively larger and 
more contorted reaction dynamics in the smaller system. 
Further reflection suggests that atoms are moving that 
would not be able to in the complete protein system, for 
example the Mg^+ ions. Wc tested this by adding two 
crystallographic water molecules to Mgi, and observing 
the reaction to become endothermic by ^ 0.5 eV. Clearly, 
the 82-atoni system is under-constrained because it is not 
large enough to provide accurate energetics of the relative 
reaction rate of the two systems. 

A potentially very useful observation is that the geome- 
tries of reaction in the 82-atom calculation for both con- 
formations share several similarities with the 244-atom 
TC calculation. This suggests that the underconstrained 
calculations can offer useful clues as to which conforma- 
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FIG. 3: (a) Energy barrier for the 82-atom system in the RC 
and TC structures. The TC transition state analog lowers the 
reaction barrier and favors the product state. The TC ener- 
gies are shifted down by 400 meV relative to the RC energies 
to facilitate comparison. Each circle represents an image in 
the NEB calculation, (b) Reaction path for phosphate trans- 
fer in RS PKA with 82 atoms. The three image sequence is 
the reactant state, transition state, and product state. The 
atoms are color-coded: Red is oxygen, light blue is carbon, 
dark blue is nitrogen, silver is hydrogen, green is magnesium 
and gold is phosphorus At the transition state the phosphate 
group is planar and the substrate proton has not yet trans- 
fered to the Aspl66 catalytic base. 



tions will favor reaction. 

An energetic offset of 400 meV has been subtracted 
from TC for ease of comparison of the barrier and 
exothermicity to RC. Although we do not place great 
confidence in this number, it in interesting to note that 
the additional interactions included in the 244-atom sys- 
tem stabilize the reactive conformation of the protein. 



DISCUSSION 

Much has been written about the reaction dynamics 
problem |23, |2J, |23, |2g, but it is primarily concerned 
with the very difficult problem of combining the vibra- 
tional transition with solvent motions. These ideas lead 
to something like |23, which docs not make use of the 
separation of energy scales which comes from understand- 
ing protein dynamics, per se, and does not provide modu- 
larity of the different aspects of the calculation (quantum 
mechanical, hydration, and conformational motions.) By 
allowing relatively few atoms to move, we increase the 
probability that the dynamics of traversing our zero- 
temperature pathway will be simple enough to occur 
quickly. The formalism of Agmon and Hopfield can then 
be used to relate the vibrational transition to an overall 
reaction rate. 

The current calculations are evaluated in the context 
of the Agmon-Hopficld model in Figure 4a, which shows 
the lATP (RC) structure at the peak of a distribution 
of conformations and 1L3R (TC) at the side. Figure^ 
shows the barrier computed in the two calculations — 60 
ksT for the reactant conformer and 6 ksT for the transi- 
tion conformer. Figure^ shows the reaction flux across 
the barrier as a function of protein conformation, simply 
the product of the probability of a conformation with the 
rate of reaction for that conformation g [cc) k (cc) , where 
k {cc) = exp [-H {cc) /ksT]. 

This model provides an estimate of the distance-scale 
over which the barrier changes can be made by com- 
bining the information in Figs. |2 and 0]:, which shows 
that a 0.9 A shift in the Mg^+ position, combined with 
several other changes of a similar size shifts the barrier 
by 50 ksT. Thus, a reasonable flux is confined to a 
multi-dimensional region only ~ 0.2 A wide! AUostery is 
the property by which small molecules or proteins bind- 
ing distant from the active site can influence activity by 
changing the protein conformation; the present calcula- 
tion shows that these changes can be quite sublte. 

The exact value of the reaction rate requires knowl- 
ege of the conformational occupancy of the low-barrier 
conformation, which we have not attempted to calculate 
in this work. There is a need for methods that can ex- 
plore conformational changes in proteins. In this regard, 
inspection of the calculation in Fig. 3 can be quite help- 
ful in discovering what conformations to look for when 
screening an MD simulation to find the correct portion 
of conformation space. 

A recent crystal structure of a Y204A mutation of PKA 
in complex with a peptide inhibitor supports two aspects 
of our calculation |28l |. First, a mixture of reactants 
and products can co-exist, and second, very few residues 
move in the course of the phosphate transfer. Yang et 
al. specifically note the absence of motion in the Mg"*""*" 
ions, coordinated water, K168, and S53 [28|. 
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FIG. 4: Schematic separation of reaction rate into confor- 
mational and vibrational (Arrhenius) coordinates, (a). The 
distribution of protein conformations, with RC representing 
an average structure and TC a less-populated member of the 
distribution, (b) We have calculated an enthalpy barrier for 
two members of the ensemble, and interpolated between in 
this figure, (c) Reactive flux as a function of conformation, 
calculated using the Agmon-Hopfield formalism, described in 
the text. 



CONCLUSIONS 

Two important lessons can be learned from this work. 
First, approximately 200 atoms need to be included to 
obtain the correct electronic structure, and also main- 
tain enough constraints on the boundary conditions to 
meaningfully reproduce the effect of protein conforma- 
tion on reaction rate. Second is the cxistance of the zero- 
temperature pathway in this particular case. In proteins 
where it is not possible to obtain TC crystal structures, 
it is essential to find appropriate methods to explore the 
variety of protein conformational motions. Since protein 
motions typically occur on timescales ranging from hun- 
dreds of nanoseconds to microseconds, it is unlikely that 
simple embedding of the quantum region in a larger clas- 
sical region will resolve this difficulty. 

Is ab-initio quantum chemistry now in a position to 
answer the biologically motivated questions posed in the 
introduction? It is clear, at least, that the active site 
cleft is not closed far enough in the RC structure. On 
the other hand, the real power of these techniques will 
only become evident when it is possible to thouroughly 
sample conformational motions of proteins and protein 
complexes. If studies of protein folding are any guide. 



this day is fast-approaching [2£ 
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